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We study steady state transport through a double quantum dot array using the equation-of- 
motion approach to the nonequilibrium Green functions formalism. This popular technique relies 
on uncontrolled approximations to obtain a closure for a hierarchy of equations, however its accuracy 
is questioned. We focus on 4 different closures, 2 of which were previously proposed in the context of 
the single quantum dot system (Anderson impurity model) and were extended to the double quantum 
dot array, and develop 2 new closures. Results for the differential conductance are compared to those 
attained by a master equation approach known to be accurate for weak system-leads couplings and 
high temperatures. While all 4 closures provide an accurate description of the Coulomb blockade 
and other transport properties in the single quantum dot case, they differ in the case of the double 
quantum dot array, where only one of the developed closures provides satisfactory results. This 
is rationalized by comparing the poles of the Green functions to the exact many-particle energy 
differences for the isolate system. Our analysis provides means to extend the equation-of-motion 
technique to more elaborate models of large bridge systems with strong electronic interactions. 



I. INTRODUCTION 

The interest in transport through conjugated 
molecules has grown in recent years in light of their 
potential applications in electronic and optoelectronic 
devices PEl While certain transport properties can be 
treated within a non-interacting picture via the tight 
binding approximation combined with the Landauer 
formalism^ often the description of transport requires 
the inclusion of many-body electron-electron and/or 
electron-phonon correlations™ For very simple and small 
systems, introducing such correlations can be done, 
for example, by means of time-dependent numerical 
renormalization-group techniques,! 5 -^ many-body wave- 
function approaches, diagrammatic techniques to real 
time pat h int egral formulationJ^J] or reduced dynamic 
methods! 12 * 13 ^ The treatment of correlations becomes 
a greater theoretical challenge in systems with many 
electronic states driven away from equilibrium] 14 ! 15 ! 
where the computational cost of numerical techniques 
increases rapidly beyond current capabilities. 

A central framework dealing with transport in large 
systems is the nonequilibrium Green functions (NEGF) 
formalism EMED The equation-of-motion (EOM) method 
is one of the more basic ways to calculate the Green 
functions (GF) of an interacting quantum system. Its 
main advantages are the simplicity and relatively mild 
scaling with the number of electrons, under simple trun- 
cation schemes! 2 ^! The EOM nonequilibrium Green func- 
tion formalism provides a qualitative description of trans- 
port phenomena in strongly correlated systems, such as 
the Coulomb blo ckade effectpiHsS] an( j ^h e Kondo effect 
in quantum dotspSEH However, questions rega rding the 
validity of the EOM approach have been raisedJ2SED p or 
example, it has been shown to violate the Friedel sum 
rule^near the Kondo regime 26 and basic Green function 
symmetry relations away from the Kondo regime.^ In 
the latter case, symmetry relations can be restored and 



at least for the Anderson impurity modelj^the approach 
recovers the Kondo peaks and provides a quantitative de- 
scription of resonant transport.^ 

In this paper we study the role of different approxi- 
mate closures to the EOM of the NEGF formalism on 
steady state properties (namely, the differential conduc- 
tance) for a double quantum dot (QD) array, coupl ed to 
two macroscopic leads (an extended Hubbard modelj^^ 
also known as the double Anderson model 32 ). Four clo- 
sures are examined; two already propose d 23 ! 33 ! (approxi- 
mation 1 and 4 described in subsection |II B 1| and |II B 4[ 
respectively) and two developed in this work (approxi- 
mation 2 and 3 described in subsections |II B 2] and |II B~3| 
respectively). The results obtained from the different 
closures were compared to the results attained using a 
many-particle Master Equation (ME) approach^ ade- 
quate for weak hybridization (system-leads couplings) 
and high temperatures) 3 ^!! I n contrast to the simplest 
case of a single site model (Anderson impurity model) 
in which different closures beyond the simplest Hartree 
approximation scheme^ yield very similar transport re- 
sults^! (steady state current and differential conductance 
as a function of the applied bias voltage) at high tempera- 
tures, we show that this is not the case for the double QD 
array, where different closures yield very different steady 
state currents and differential conductance curves. 

The performance of the different closures is analyzed 
in terms of the poles of the GFs in comparison to the 
exact many-body result for the isolated system. We find 
that one of the closures developed in this work provides 
the most accurate description of the poles and also the 
best overall agreement with the ME approach for all pa- 
rameters studied in this work. While these results are 
encouraging, a word of caution is in place. It is clear 
that the conclusions drawn from the performance of the 
different closures for the single site model cannot be ex- 
tended directly to the two site model, in analogy, a suit- 
able closure for the two site model may fail in the larger 
systems. Thus, the study of larger arrays of QDs will 
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Figure 1: A sketch of the double QD bridge. See main text 
for the definition of all quantities. 



require analysis along the lines sketched in this work. 

The paper is organized as follows: in Sec. [IT] we present 
the double QD model Hamiltonian, provide a short de- 
scription of the equation-of-motion technique and a de- 
tailed description of the different approximate closures 
to the EOM. Sequential labeling of the different approx- 
imate closures refers to the order of the closure. Results 
and discussion are given in Sec. |III| for the cases of the 
symmetric and asymmetric bridges. In Sec. |IV| we con- 
clude. 



II. THEORY 

A. Model Hamiltonian 

We consider a system of two coupled QDs array con- 
nected to two macroscopic leads, as sketched in figure [T] 
The Hamiltonian has the following general form: 



H — Hb + Ha + Hi 



(1) 



with Hb describing the macroscopic leads (left and right 
contacts), H$ describes the system of interest, and Hi 
is the interaction Hamiltonian between the system and 
the leads. The leads (left (L) and right (R)) are modeled 
as infinite non-interacting fermionic bathspSEU anc j are 
assumed to be each at its own equilibrium, characterized 
by chemical potentials fi^ and /i^, where the difference 
Mi — M-R = e ^ is the applied voltage bias. The leads' 
Hamiltonian is given by: 



H, 



a,k£{L,R} 



■ Cko 



(2) 



where 6ka is the energy of a free electron in the left or 
right lead, in momentum state fc and spin a. The oper- 
ators Cfccr/cl are the annihilation/creation operators of 
such an electron. The double QD system is described by 



an extended Hubbard modelf^EH 

Hs = ^ s m<J n m(J + U n m fn m i 

+V ^ n aa np a ' +h^2 {dtadpa + h.c.) , (3) 

a, a' o 

where n aa — d a(7 d acr is the number operator of the elec- 
tron occupying site (dot) a with spin a and energy e aa , 
U is the repulsion energy between two electrons on the 
same site with opposite spins (intra-dot repulsion), V is 
the repulsion energy between two electrons on different 
sites (inter-dot repulsion), and h is the coupling strength 
for electron hopping between the two sites. The interac- 
tion between the system and the contacts is simply given 
by the tunneling HamiltonianPSl 

Hi=J2 (*te C L < k»+' l - C -)+ £ (tlAadpa + h.C. 



a,k£L 



cr.keR 



(4) 

The parameter i£ m represents the coupling strength (hy- 
bridization) between the system and the leads, and the 
index m runs over the site index {a,/3}. 



B. Equation of motion 

The above model consists of an interacting system cou- 
pled to two electron reservoirs with specified chemical 
potentials and temperatures. In order to obtain a solu- 
tion to this many-body out of equilibrium problem, we 
res ort to t he EOM approach within the NEGF formal- 
ism! 4 ! 25 ! ^ We begin by defining the contour ordered GF, 



, x T c ^ a (t 2 ) % (ti)), where tf/tft are 



the system's annihilation/creation field operators, and 
Tc is the contour time ordermgoperatorPESI The EOM 
for the contour ordered GF 43 * 44 ! is obtained from the 
Heisenberg EOM for a Heisenberg operator ^-A(t) = 



H(t),A H (t) 



A, B 



= AB-BA. 



f t A H (t) , where 
A full description of the system requires the knowledge of 
the retarded, advanced and lesser (distribution) GFsP^ 



Gap (h,h) 
G a p {t2,h) 

G af3 (h,tl) 



0(ta-*i)({$a (*a)»*J (*i)} 



(5) 



where {A, B} is the anti-commutator of A and B. These 
real time GFs can be extract ed fr om the contour ordered 
GF using the Langreth rulesP^l Except for very simple 
models (e.g. see Refs.E7|and|I5) , the EOMs of the NEGF 
will produce "new" and higher order GFs that need to be 
evaluated. In general this leads (after a few iterations) 
to a non tractable hierarchy of equations. To obtain a 
working closure one has to truncate the resulting set of 
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equations and/or decouple the higher order GFs and ex- 
press them via lower order ones. Recently it has been 
shown that such a procedure may lead to symmetries vi- 
olations that the GFs must obey by definition^ and a 
routine to restore back the symmetries was suggested.^ 
In what follows we examine 4 different closures to the 
EOMs of the NEGF for the Hamiltonian discussed in sub- 
section |II A[ and apply the symmetry restoration scheme 
to circumvent the inherent flaw of the EOMs approach. 
Our starting point is the contour ordered GF: 

GZ ( f,(t,t') = -l(T c d aa (t)dl (r (t')). (6) 

Writing down its EOM and expressing the result in 
Fourier space (as we are interested in steady state prop- 
erties) we find that the single particle GFs obey: 

x (6 afj + hGpp H + VG%% M 
+VG%£ M + UG%£ («)) , (7) 

where E° CT (w) = J2k l^fcJ 2 ~ £ fco-) _1 is the tunneling 
self energy (the label "0" refers to the limits U, V — > 
where the total self-energy is given by TP aa (w)) and 
^afhf ( w ) ^ S t ne Fourier transform of the 2-particle GF 

= -i{Tcn aT {t)dp a (t)d\ a {t'))- Calculat- 
ing the EOMs for the higher order (2-particle) GFs will 
give rise to other 2-particle and 3-particle GFs: 



and 

x (n ai<T ) - 5 ai (<fl a<r dp t<T )) 
+hGJ™ ( W ) + UG™% (w) 

+ E(«^7 H-CF^ (w)) 
k 

+ E**7 F ^7 M + ^ G ^7 M ] -( 10 ) 

In the above equations, G^f (w), (w), (w), 

F afc/37 ( w ) an( ^ F fea^7 ( w ) are ^ ne Fourier transforms of 

Gj^(M') = -^(2b^(*)^(«)i r «r(*)4 ff (* / )) ) 
(*, O = -| (T c n ar (t) n Ps (t) d ia (i) d\ a (f )), 
F^ 7 (M') = -|(T c n QT (t)c feCT (t)4 ff (f)>, 

TO (*' = "I < T ^- (*) (t) d^ty (t) d\ a (f )) 
and F££ (i, f ) = -i (Tc4 (i) d QT (f) (t) dt CT (f )) 
respectively. To continue, one has to formulate the 
equations for these new GFs, which in turn will lead 
to other (higher order) GFs. This infinite hierarchy of 
equations needs to be truncated at a certain level, a 
process which is referred to as "closure". In general, clo- 
sures cannot be improved systematically. Furthermore, 
it is often difficult to assess, a priori, the accuracy of a 
given closure. We now discuss several different closures 
which are physically motivated, tractable, and some are 
commonly used in the context of transport. 



(nu- Eptr -v-j^ a H) _1 

x (S yf} {n ara ) + hG%% H 
-hG$™£ H + ftG*£ ( W ) 

+ E(M!^ H-CF^ 7 (w)) 
fc 



A- 



;^F^ (w) + C/G^ (w) 



7m M) , 



(8) 



(/^ - - (7 - E° CT (w))" 1 

X (<5 7Q (V> + ^^7 M 

-M3^ («) + M3»^ (w) 



-E 



-7G^7 M) . 



(9) 



i. Approximation 1 

Following the derivation given in Ref. 2.'? the fol- 
lowing approximations are made: (a) all 3-particle 
GFs are set to zero, (b) simultaneous tunneling 
of electrons of opposite spins are neglected, (c) 
GFs mixing leads and system operators are decou- 

pled scPF(M') = -i(r c c ka {t)n aT {t)dl a {t')) « 

-Rt* 7 /d*lff*(*.*l)( I bd 7 a(*l)nar(tl)4«r(* / )) where 

(iftjj — £fccr) fffe (Mi) = S (t — ti), which is obtained by 
assuming that n aT (t) is constant (as is the case in 
steady state). (d) The remaining 2-particle GFs of 
the form G££*(*,i') are decoupled so G^(t,f) = 

(n ar (t)) (t,t'). Assumption (c) is equivalent to 
treating the coupling to the leads up to the second or- 
der with respect to tZ _. It neglects processes necessary 
to qualitatively capture the Kondo effect J^H^ yet results 
are predicted to be reliable for temperatures above the 
Kondo temperature (Tk The resulting equations 
are given by (for brevity, we omit the implicit depen- 
dence on ui): 
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G° 



G 



Pa 



X(l 

(hw - ep a 
x (hGZ - 



•4 <J(J(T 

u aaa 



VGfcl + VG$Z) > 



(11) 



(12) 



In general the GFs depend on the expectation values of 
(n 7T ) = (d^dyT^ and (d^dp^^, which are determined 
by means of the lesser GF: 

■ J- />00 

( d t^>=-^ y jg^m)^, (is) 

thus, a self consistent calculation is required. 



,8. Approximation 2 

A seemingly better approximation scheme is one that 
relaxes the last mean-field approximation (assumption 



"d") described in subsection II B 1 and the 2-particle GFs 
of the form &™£ (t,f) = -| <T c n QT (i) (t) < CT (O) 
are not decoupled but treated fully. Thus, while equation 



(111 remains 


the 


same, the equations for the 2-particle 


GFs will now 


be j 


jiven by: 




(h0J - Spa - 


V - 


~ ^P<7> ^aPa 


= hG^aa, 


(hw - ep a - 


V - 


V>0 \ (T^ (7(7(7 

^Pa) ^aPa 


= h<G Paa - ( d la d P,cr) , 


(hw - ep a - 


-Ti- 


vO \ (P (7(7(7 
~ ^Pa) ^PPa 


= hGpaai 


(fruj - e aa - 


ll - 


V-iO \ /f* (7 (7 <7 

~ ^aa) "aao 


= ("a?) + hG° a pa, 


(hw - £ aa - 


V- 


VO \ (T*(7(7(7 


= K) + /.G^, 


(hw - £ aa - 


V- 


VO \ /P (7(7(7 

^aa) ^Paa 


= (n^) + WJ|^. (14) 



As noted above, the retarded, advanced and lesser GFs 
can now be evaluated using Langreth theorenPS an£ } the 
expectation values ((n 1T ) — ^<i 7T G? 77 -^) and (d^dp^) are 
determined via the lesser GF. 



3. Approximation 3 

A more complete treatment of the 2 nc ^ order GFs re- 
quires relaxing assumption "b" in addition to assump- 
tion "d" (see previous subsection), described in subsec- 
tion |II B 1| Namely, we attend to the 2-particle GFs 



that describe simultaneous tunneling of electrons of op- 
posite spins in the double QD system, G^p^p (^ = 

-i (r c di s (t) dps (t) d aa (t) d\ a (i')). Again, the only 

changes are in the equations for the 2-particle GFs, and 
the resulting equations are given by: 



1 (7(7(7 

"apa 



V 



(hw 


— Spa- 


-u - 




lP (7(7 (7 

^PPa 


= h(np s )G°° a -, 


rr^GGG 


(hw 


-£p* 


— V — 




/p (7(7(7 

^"aPa 


= /i(n^)G^-(4^ CT > ) 


^a/3a 


(hw 


— £pa 


-V- 




fP GOG 

^aj3a 






(hui 


£-Ct(7 


-u- 




fP (7(7(7 


= (naff) + ft (flaa) G^, 


fP ctctct 

^,3 /3a 


(hw 


£-aa 


— V — 




/p (7(7(7 

^"paa 


= (npa) + h (n aa ) Gp°, 




(hw 


£acr 


— V — 






= (np s )+h(np s )G°p° a . 


/p <TfT(X 



L(r<(7(7(7(7 

n^uppa 



= {huj - sp a 
= (tou-epr-U-TPfr) 



hG$l%) 



x WE?SS 



n^paPa 



hG^pZ) 



-U-Y. 

hG%l 



h(T aac " 7 



1 (717(7 

"Paa 



i (7(7(7 

"Paa 



T /3a/3a 



ff /3aaa 



v aP pa 



u aPaa 



X {(naa) ~ 
(Hu-Ear-V 

x ((np a ) + hGl%) , 

x ((np s ) + hGfpl + hG 
(hw - ep a - e aS + ep s ) 
x (hG%l - hG%l + hG% 
(hw — e aa — e aS + ep s ) 

x ((4^)-hGZl- 
(hw - ep a - sp s + e aS )~ 



Paaa 
1 



hG 



a/3aa ) 



a. a) i 



hG$Z 



j^/p era" (jo" 
n ^PaPa ) ' 



(7(7(7 U(T*(7(7(7 

aPa ~ n ^ppa 



•i G(7(7(7 \ 

11 aPaa j 



— (hw — E aa — Ep s + Eaa-) 
X {(d^dpa 



hGZa 



■hG 



(7(7(7 

aaa 



hG^pfa 



(15) 



This case is not different from the previous two in the 
sense that a self consistent treatment is required. 



4- Approximation 4 

Finally we follow the derivation of Ref. 33; Here, the 
following approximations are made: (a) simultaneous 
tunneling of electrons of opposite spins are neglected, 
(b) GFs mixing leads and system operators are de- 
coupled so ¥(t,t') = -i(T c c k< ,(t)n aT {t)d\ a (t)) w 

-it^I^hgkit^iTcdjait^riari^d^if)) (see 

discussion in subsection II B 1| , and (c) 3-particle GFs 
are introduced in a mean-field way, i.e. 3-particle GFs 

of the form -\ (r c n lcr ( (t) n$ T (t) d aa (t) d\ a (t')\ are 
decoupled to -| (n 7CT / (t)) [T c ns T (t) d arJ (t) d^ (t')\ - 
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I {n Sa (t)) ^T c n 7CT , (t) d QCT (i) 4, (*') 
equations are given by: 



The resulting 



^ CTCTCJ" 



•< (J<7<J 



i (7(7 CT 



•> (T(T<7 

u aaa 



Hid 



h&%™ + {np s ) V 



a 

pact 

hcu — S aa — U 
{np&) + hGffi a + {nps 
tvjJ - e aa - U (n aS ) - V {{nps) - 
M + h&™p° a + {np„) (UG™ 
hoj ~ ep a ~U {nps) - V ((n„) - 
+ (n^) (UG 

hu - e aa - U - V {{nps) + {np a )) - I 
{n a s) + hG%l + {n a s)V{GpZl + G% 

hoj -Bpcr-U {nps) - V {{n aS ) + 1) - 



4.) 



U ~V{{n aa ) + {n a s)) 

9 )-V{{np <T )+l)-Y,l a )-' 

■VG$Z)), 



1) 



VG a a %)) , 



y{) 



r 
)). 




Bias Voltage (O/U) 



Figure 2: (color online) Plots of the differential conductance 
versus the bias voltage for the symmetric bridge ( e a f = £ a \ = 
£/St = e /3J. ~ 0.35(7) for V = 0. Upper left, upper right, lower 
left and lower right panels correspond to h — 0.1(7, 0.3(7, 



As before, the different expectation values {{n 1T ) = 
(d 7r d 7T ) and (d^dp^y) are obtained by integrating the 
relevant lesser GF, which is calculated from the contour 
ordered GF. 



III. RESULTS AND DISCUSSION 



VGZpa)) -°- 5U and °- 7C/ > respectively. Black curves corresponds to re- 
/i^sults based on the ME. Red (circles), green (diamonds), blue 
(triangles) and magenta (stars) correspond to the results ob- 

Thc 



tained by approximation schemes 1 to 4, respectively, 
notation |z) — > \j) indicates that the conductance peak calcu- 
lated by means of ME corresponds to a transition form the rij- 
particle states to any of the n^-particle states. The remaining 
model parameters were V^ La = F^ a = F^ = F^ = 0.015(7, 
t _ r4 



= = rL = r* Q = o, and /r 1 = u/4o 



-4 

L0 



At this point we wish to examine the different approx- 
imations and find which one leads to a system GF that 
describes the double QD more accurately or at least qual- 
itatively. As we are interested in transport properties of 
a system weakly coupled to the macroscopic leads we will 
compare the results obtained from the different approx- 
imations to the ones calculated using the many-particle 
ME approach.^ Under th ese ass umptions, the ME is be- 
lieved to be fairly accurate. ^ 53 ! We wish to note that the 
EOM approach is not limited to the weak coupling case. 
As a measure of the quality of the approximations we 
chose to calculate the differential conductance, dl/d&. In 
the many-particle picture, in the wide band limit (where 
the interaction with the leads only broadens the energy 
levels of the system without introducing any spectral 
shift), we expect to observe peaks in the differential con- 
ductance at values of the bias voltage that correspond 
to n L/R -E f = ±e| « AE (N) =E(N)-E(N- 1), 
where Ef is the equilibrium Fermi energy of the elec- 
trodes (throughout taken to be zero) , E (TV) is the en- 
ergy of the many-particle state with N electrons of the 
unperturbed system, and $ is the applied voltage. The 
voltage can be applied symmetrically to both leads (i.e., 
fiL = Ef + and /ir = Ef — e%), or asymmetrically 
{fiL = Ef + e& and = Ef). In the present study we 
have used the symmetric version. 

The differential conductance is derived from differen- 
tiating the steady state current with respect to the bias 



voltage and was evaluated from the Meir-Wingreen for- 
mula!^ 

1 = £hj d£ (^ifL^-^)r L {e) 

x (G r (e)-G a (e))}+Tr{T L G<(e)}).{17) 

In the above, G r (e), G a (e) and G K (e) are the the re- 
tarded, advanced and lesser GFs, respectively and Yl is 
the matrix coupling of the system to the left reservoir, 
with elements (r L )* Q = 2tt £ S (s - e ka ) |i£j 2 . The 

resulting EOMs were solved self-consistently in Fourier 
space with a frequency discretization of dcj = 0.0005(7 
over 32, 768 grid points. Depending on the approxima- 
tion, 15 — 150 self-consistent iterations were required to 
converge the results. Convergence was declared when the 
population values ((n mT )) at subsequent iterations did 
not change within a predefined tolerance value chosen to 
be 10 -6 . For each set of calculations symmetrization rou- 
tine was applied to restore the symmetry relations of the 
GFsP 



A. Symmetric bridge 

The transport through the double QD system can be 
classified into symmetric and asymmetric bridge setups, 
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Bias Voltage (O/U) 

Figure 3: (color online) Same as figure [5] but for V = Q.W. 



GF poles 


Energy differences 




AE(1) : 


e — \h\ 




e + \h\ 


e±|ft| 




AE{2) : 


e + ^0 - Si 




e+(U + V) /2-S 2 ±\h\ 


e + V-\h\ 


e+(U + V) /2 + S 2 ±\h\ 


e + V+\h\ 


e + V ±\h\ 




e + U±\h\ 



Table I: Left column: Location of the poles of the unper- 
turbed system's GF as calculated using the 2 nc ' approxi- 
mation. Right column: The differences in energy between 
many-particle states that differ by one electron, such that 

AE (N) =E(N)-E(N-1). HerSi = § \J {U - V) 2 + Ah 2 , 
and S 2 = \ \]{U - V) 2 + 16/i 2 



with or without inter-dot repulsion term, V. In this sub- 
section we first consider the symmetric setup in which 
e a -f = E a i = s ^ = Epi = e = 0.35U. The remain- 
ing model parameters were taken to be rt = ri = 

r T ^ = = 0.015C/, = = rL = rj^ = o, 

and = [7/40. In figure [2] we plot the differential 
conductance for a symmetric bridge for different values 
of the hopping term h. We set the inter-dot repulsion 
V = 0. The black curves (solid line) are the results ob- 
tained by fully diagonalizing the bare system H s (and 
solving the ME) . The other curves represent the outcome 
of the NEGF formalism within the different closure ap- 
proximations. We also label the different peaks in the 
differential conductance, obtained via the ME approach, 
with the corresponding transitions between many-body 
states, i.e., |0) — > |1) corresponds to transitions from an 
empty system to a system with a single electron, etc. In 
the single particle GF formalism, peaks in the differen- 
tial conductance will occur at the poles of the calculated 
GF. Hence a good approximation is one that will pro- 
duce single particle GF with poles at the position of the 
many-particle transitions. 

For the smallest value of h, we find that all approxima- 
tions agree qualitatively with the ME approach. When 
the value of h is increased it is clear that approximation 
1 (red circles) breaks down, implying that this simple 
closure is insufficient to describe strong hopping between 
the quantum dots. Approximations 2 and 3 (green di- 
amonds and blue triangles, respectively) do agree with 
the ME, but "miss" certain conductance peaks (e.g. the 
peaks at */[/ pa 0.75 in the upper right panel, */[/ ps 0.5 
in the lower left panel and the peak at */[/ ss 0.2 in the 
lower right panel) , all of which correspond to transitions 
involving a 2-electron occupancy. Approximation 4 (ma- 
genta stars), which includes 3-particlc GFs at a mean- 
field limit, performs slightly better in this respect. 

In figure [3] we present results for the differential con- 
ductance obtained for the symmetric bridge where the 



inter-dot repulsion, V, is included. All other parameter 
are identical to those of figure[2] In this case, we find that 
approximation 1 is not suitable even for small values of 
the hopping term h while approximation 4 appears to 
work for low values of h < | U (both upper panels) , how- 
ever, it fails to capture peaks resulting from transitions 
through 2-electron occupancy at higher values of h, as 
depicted in the lower panels of figure [3) We note that for 
this set of parameters, such transitions involving 2 elec- 
trons are absent for h < |£7 in the bias voltage studied. 
Approximations 2 and 3 agree very well with the ME re- 
sults for all values of h, even at values of the bias voltage 
that correspond to transfer through 2-electron states, in 
contrast to the case where V = in which they fail to 
capture conductance peaks involving 2 electrons. 

The performance of the different approximations can 
be rationalized in terms of the pole structure of the un- 
perturbed system GF, which can be compared to the ex- 
act many-body energy differences between many-particle 
states that differ in one electron. While it is possible to 
carry out this analysis for all closure approximations, it 
is often a tedious task. Thus, in what follows we pro- 
vide such analysis for the case of approximation 2 only. 
The poles of the GF, the many-particle energy levels and 
the differences in energy are summarized in table [T] (see 
Appendix for more details regarding the derivation of 
the poles of the GF). As can be seen from table [I] con- 
ductance peaks corresponding to transitions |0) — > |1) 
and |1) — > |0), that is, peaks appearing at the values 
of AE(1), are captured by approximation 2 since the 
GF has poles at the correct locations. Higher excita- 
tions involving 2 or more electron occupancies are not 
fully or systematically accounted for by approximation 2 
(or any of the other closures described in this paper, for 
this matter). In general we find that such higher transi- 
tions are not captured by approximation 2 when V <C U . 
For V = U we find that AE (2) has 4 different values: 
e + V ± \h\ and e + V ± 2 \h\, concurrently the GF has 
poles at e + V± |h|, thus, some of the transitions involv- 
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ing the 2-electron states \N — 2) (particularly those with 
E(N = 2) = e + V± \h\) are described by the NEGF. 

Following this short analysis we can now better ex- 
plain the results of figure [3] For a large value of the 
inter-dot repulsion (V = 0.8t7, thus V ~ U), one ex- 
pects that the calculated GF will capture the higher or- 
der transitions in the relevant bias window and agree 
with the ME results. If one considers the second peak 
rj 1.1) in the lower right panel of figure [3] it re- 
sults from transmission through a many-particle level 
with \N — 2). For the symmetric bridge, this peak cor- 
responds to AE = (2s +§(£/ + V) - S 2 ) - (e - \h\) = 

e+UU + V)-S2+\h\, where S 2 



^(u-vy + wh 2 . 

Under the assumption that V ~ U, one can approxi- 
mate AE m e + V — \h\, whereas the GF has one of 
its poles at Pq — e + | (U + V) — Sx, where Si — 
l 

2 



-V) 2 + 4h?, 



which, for V ~ U can be approx- 
imated by Pq ps e + V — \ h\. For the studied parameters 
(see figure [l| we find P G = 0.54289 and AE = 0.54643, 
and indeed the differential conductance based on approxi- 
mation 2 show a peak at twice this value*/ u as 1.1. While 
for the case where V = this transition is overlooked. It 
is easy to verify that a similar argument holds for the 
second peak in the lower left panel of figure [3] as well. 



B. Asymmetric bridge 




1 1.5 0.5 
Bias Voltage (O/U) 



1.5 



Figure 4: (color online) Plots of the differential conductance 
versus the bias voltage for the asymmetric bridge {e a f = 
e al = 0.15U and s Pt = e pi = -0.2(7) for V = 0. Upper 
left, upper right, lower left and lower right panels correspond 
to h — 0.1U, 0.3U, 0.5U and 0.7U, respectively. Black curves 
corresponds to results based on the ME. Red (circles), green 
(diamonds), blue (triangles) and magenta (stars) correspond 
to the results obtained by approximation schemes 1 to 4, re- 
spectively. The notation \i) — > \j) indicates that the con- 
ductance peak calculated by means of ME corresponds to a 
transition form the n^-particle states to any of the n^-particle 
states. The remaining model parameters were F 



t 



1 R0 



1 Rf> 



0.015U, 



/J- 1 = U/40. 



L L/3 



1 La — 

0, and 



We now turn to discuss the case where e aa ^ ep a re- 
ferred to as the asymmetric bridge. Once again we have 
calculated the differential conductance using the 4 differ- 
ent closure approximations to the NEGF formalism and 
compared the results to the differential conductance ob- 
tained by the ME. Analysis based on analytic expressions 
for the poles of the GF or the many-particle energies of 
Hg is more difficult, and the expressions are not as com- 
pact as in the symmetric case. The results for the poles 
of the GF within closure approximation 2 are given in the 
Appendix, while the many-body energy differences were 
obtained numerically. 

In figures [4] and [5] we plot the differential conductance 
for the asymmetric bridge for different values of the hop- 
ping term h for V — and V — 0.8U, respectively. The 
on-site single particle energies were e a ^ — e a ^ — 0.15J7, 
£pt = £ Pi = — 0.2[/ . The remaining model parame- 
ters are identical to those of the symmetric bridge and 

= = o.oisu, 

0, and p- 1 = U/40. As be- 
fore, the black curves (solid line) corresponds to the ME 
results. The other curves represent the outcome of the 
NEGF formalism within the different closure approxima- 
tions. We also label the different peaks in the differential 
conductance with the corresponding transitions between 
many-body states, i.e., |0) — > |1) corresponds to transi- 
tions from an empty system to a system with a single 
electron, etc. 



were taken to be 1^ 



1 Lc 



■p4- 

i L/3 



1 Ra 



1 Ra 



From figures [4] and [5] it is obvious that approxima- 
tions 1 and 4 do not perform as well as approximations 
2 and 3. We would like to note that approximations 1 
and 4 utilized a mean-field like approximation decoupling 
the higher order GFs, while in approximations 2 and 3 
higher order GFs are ignored altogether. For all parame- 
ters studied in this work (not all presented here) , we find 
that approximation 2 performed better than all the other 
approximations, suggesting that including higher order 
correlations in a mean field fashion or a more complete 
treatment of the 2 nc ^ order GFs is not advantageous. 

We find that approximations 2 and 3 predict negative 
differential conductance at higher values of h. not ob- 
tained by the ME, as shown in the lower panels of figure 
|U The dips occur (in both cases) at values corresponding 
to the activation of the anti-bonding single electron state. 
While this transition is suppressed in the ME approach , 
it appears to be enhanced in the NEGF formalism. 



IV. CONCLUDING REMARKS 

In this work we have assessed the validity of the EOM 
approach to the NEGF formalism for an interacting sys- 
tem coupled to two macroscopic leads. The interacting 
system consisted of two coupled quantum dots, each with 
one electronic level (spin up and spin down), connected 
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(Anderson impurity model), they are not easily expand- 
able to systems with more complex noncquilibrium dy- 
namics, such as the double QD model (double Anderson 
model). In light of this, the success of approximation 2 
for the double QD model does not necessarily imply that 
it will provide quantitative results for a more elaborate 
system. However, analysis of the poles of the resulting 
bare GF and comparing them to the exact many-body 
energy differences does provide a tool to assess the accu- 
racy of a given approximate closure and can be used for 
larger bridge systems even if these can only be performed 
numerically. 



0.5 1 1.5 0.5 1 1.5 
Bias Voltage (4>/U) 



Figure 5: (color online) Same as figure [I] but for V = 0.8(7. 
Results obtained from approximation 1 are only presented for 

the case h = Q.1U (upper left panel) as we could not converge V. ACKNOWLEDGMENTS 

it for higher values of h. 



serially, taking into account intra and inter-dot Coulomb 
interactions. 4 different closure approximations to the 
EOM, some are commonly used in the literature and oth- 
ers have been developed here, were examined. As a mea- 
sure of the quality of the approximations we calculated 
the differential conductance (derived by differentiating 
the steady state current with respect to the bias volt- 
age) and compared the results to those obtained by the 
ME approach, which under the approximations of weak 
coupling to the leads and high temperature provides ac- 
curate results. Two different cases corresponding to a 
symmetric bridge (e aa — sp a ) and an asymmetric bridge 
(^aer 7^ £/3cr) with and without inter-dot Coulomb repul- 
sion (V), were studied for different values of the inter-dot 
hopping term h. As expected, we find that keeping more 
terms in the closure or including higher order correlations 
in the EOM, does not necessarily improve the approxi- 
mation. As a rule of thumb, neglecting higher order GFs 
(approximations 2 and 3) preforms better as compared 
to closures that include such terms at a mean-field level 
(approximation 1 and 4). 

To assess the performance of the different closure ap- 
proximations, we compared the pole structure of the un- 
coupled GF with the exact results of the many-particle 
states. Focusing on approximation 2, which provides 
the overall best agreement in comparison to the ME ap- 
proach, we found that transitions involving only single 
electron states were reproduced by the NEGF. However, 
when higher many-particle states are involved, the accu- 
racy of the approximation depends on the strength of the 
Coulomb coupling V. In cases where V ~ U, approxima- 
tion 2 also captures conductance peaks associated with 
transitions involving two electron states. 

While all the approximations described in this work 
capture the Coulomb blockade and the main characteris- 
tics of the steady state transport for the single QD model 
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Appendix 



Using the assumptions described in Sec. |II B 2[ the 
EOMs for the NEGF of the unperturbed system are (in 
Fourier space): 



G 



OCT 

a a 



0a 



(hu - e aa ) 

x (1 + hGZ + UGZ a a + VG%Z + VG$Z) , 
(hu - ep a )~ l 

x (hGZ + UGffe + VGZZ + VG%1) , 

(18) 



(hu 


— £/3cr 




— hGZZi 




(hu 


— £/3cr 


-V)G™f a 


= hGZl- 


- (dLdp 


(hu 


— £/3cr 


-U)G$Z 


= hGZ a a i 




(hu 




-U)GZl 


= (n a s) + 




(hu 




-V)GZ a a 


= ("/3a) + 


hGZZ i 


(hu 




-V)GfZ 


= (npa) + 


hGpZ- 



(19) 
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It is clear that the equations for the 2-particle GF close 
among themselves, so a simple substitution yields: 



Gcraa 
a/3a 



Gcrtra 



G" 



aaa 



Gcracr 
/3aa 



G (7(7(7 
0aa 



Hw — ep a -V ■ 
h (n aS ) 

(Hw - E aa - U) 

Hw - ep a - v - 



(Hw - e aa - U) 



h 2 

(Hw - e aa - V) 



Hw — Ep a — U 



(Hw - e aa - V) 



h (np&) 



(Hw - e aa - V) ' 
Hw — e aa — U 
x (n a a) 
Hw - E aa -V 



h 2 



(Hw - Ep a - V) 
h 2 



(Hw - ep a - v) 



Hw — e nrT — V — 



(Hw - Ep a - V) 
h 2 



(Hw - Ep a - U) 



G (7(7(7 
apa 



G (7(7(7 
a/3a 



h (n aS ) 



G a 



Go 



i a a a 

*paa 



G (7(7(7 
0aa 







h 2 ) 




a j (nw - 


T \ 

Jjcxv J 




^QU ) 


- n ) 


n (Ho- 


\ 
/ 




/ / hy. i — a- \ { nil) — 


- X )- 


h 2 ) ' 


In -) (fin} - 
\ ll a>(T/ \' 1>JJ 


-To \ 




( ( hy. ) — 7* || h/ , j — 


- Xo ) - 


h 2 ) 


(n^) (Hw - xp v ) 




■d/3,a) 


((Hw — x av ) (Hw 


- xpv) ~ 


-h 2 ) 


(np s ) (Hw - 


" X/3 U ) 




((HUJ — X av ) (Hw - 


" Xff u ) - 


h 2 )' 



(23) 



We now substitute the set of equations ( 23 1 into equa- 



tions (22) 



(20) Finally 



GZ = \(hu - x a ) - 



(hw — xp) 



x 1 + 



(hw — xp) (huj — xp) 



(hxl — Xp) 

i a a a 

"Paa 



"paa) 



(24) 



Define: 



C a/P 
v I Pv 
a/ Pu 



and rewrite the EOMs 



G 



OCT 
MM 



Pa 



(fajJ — X a ) 1 

x (1 + hGf a 
(Huj — xp) 1 
x (HG\ 



^aa j Pa ) 
£ ao/p<7 + V, 
£ n(7/Pcr + U, 



UG, 



<JOO 



(21) 



VG%£ + VG$Z) , 



-1(7(7 

1 aa 



VG^pl 



VG%1) 



(22) 



a a 
a a 



((Huj — x a ) (Huj — xp) — h 2 ) 
h 2 U(np Sl 



x I I + 

+ 



((Hut - xp u ) (Huj - x av ) - h 2 ) 
h 2 V(np a ) 



((Hw - xp v ) (Huj - x av ) 


-h 2 ) 


HV (S aa dp^) (Huj-x 


av) 


((Huj - xp v ) (Huj - x av ) 


-h 2 ) 


h 2 V (n aS ) 




((Huj ~ xp v ) (Huj - x au ) 


-h 2 ) 


U (n aS ) (Huj — xp v ) (Huj 


-xp) 


((Huj - x au ) (Huj - xp v ) 


-h 2 ) 


(npa) V (Huj - xp) (Huj - 


- Xp v ) 


((Huj - x av ) (Huj - Xp v ) 


- h 2 ) 


hV(dl, a dp <a ) (Huj - a^) 


((Hw - x av ) (hu - xp v ) 


-h 2 ) 


(npa) V (Huj — xp) (Hw - 


- xp u ) 


((Hw - x av ) (Hw - xp u ) 


-h 2 ) 



(25) 
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From the last equation we see that the NEGF has poles 
at 



(hui 


— x a ) (hu> - 




-h 2 


= o, 


(huj — 


x av ) (huj - 


X/3v) 


-h 2 


= o, 


(hw- 


x au ) (fkj - 




-h 2 


= o, 


(hw- 


x av ) (hw - 


Xf3u) 


-h 2 


= o, 



(26) 



or equivalently 




pi, 2 1 

~ 2 


(x a + x fj )±\J (x a - x p ) 2 + 4h 2 


i 


P 3,4 _ 1 

2 


(x av + x fjv ) ± \J (x av - x Pv ) 2 + 4/i 2 


p5,6 1 

~ 2 


(x av + x,3 U ) ± \J (x av - x,3 U f + 4/i 2 


p7,8 1 

2 




Ah? 



(27) 
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